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ABSTRACT 

We investigate a simple model for globular cluster (GC) formation. We simu- 
late the violent relaxation of initially homogeneous isothermal stellar spheres and 
show that it leads to the formation of clusters with radial density profiles that 
match the observed profiles of GCs. The best match is achieved for dynamically 
unevolved clusters. In this model, all the observed correlations between global 
GC parameters are accurately reproduced if one assumes that all the clusters 
initially had the same value of the stellar density and the velocity dispersion. 
This suggests that the gas which formed GCs had the same values of density and 
temperature throughout the universe. 

Subject headings: globular clusters: general — methods: A^-body simulations 



1. INTRODUCTION AND MOTIVATION 

Globular clusters (GCs) are massive, dense clusters of stars, and in many galaxies they 
are the oldest datable objects. They contain important clues to early star formation in the 
universe, and to the formation history of our Milky Way and other galaxies. The dynamical 
evolution of GCs has been studied for many decades, and is considered to be well-understood. 
However, the formation of GCs is still an open question. 

The internal structure and dynamics of Galactic GCs are well described by a family of 
lowered isothermal models (so-called King models. King 1966). The radial profiles of King 
models have a flat core (characterized by a core radius tq) and a sharp cut-off at a tidal 
radius rt. GCs are presently in dynamical equilibrium. Initially they could have formed in 
some non-equilibrium configuration which relaxed to a King model through violent relaxation 
(Meylan & Heggie 1997) followed by a slow evolution in the galactic tidal field. 
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Currently no models of GC formation explain all the observed properties of these star 
clusters (Djorgovski & Meylan 1994; McLaughlin 2000). Correlations between parameters 
which are not significantly affected by the dynamical evolution of a cluster (such as half- 
mass radius, central velocity dispersion and binding energy) are thought to reflect the initial 
conditions for GC formation. Therefore, any model of cluster formation should aim to 
address the observed correlations. 

In this paper, we investigate the simplest possible model for GC formation: the dy- 
namical evolution of a homogeneous isothermal stellar sphere. Simulations of the collapse of 
homogeneous stellar sphere with constant velocity dispersion were performed twenty years 
ago (van Albada 1982) with the intention of reproducing the r^^^ density distribution of 
elliptical galaxies. Those efforts were not very successful — the collapse produced instead a 
core-halo structure, somewhat like a GC. In this paper, we compare a family of such models 
with the observed density profiles of GCs, and attempt to use the family of such models to 
reproduce correlations between parameters of GCs in the Milky Way. 



2. MODEL 

We propose a homogeneous isothermal stellar sphere as an initial non-equilibrium con- 
figuration for GCs. To test this idea, we evolve stellar spheres with a total mass M, an 
initial density pini and an initial velocity dispersion (Jini for tend = 10 — 5000 initial crossing 
times tcross = {Rfni/^Y^'^ 1 where i?ini = (3M/47rpini)^^^ is the initial radius of the system. 
(In this Letter we assume a system of units where the gravitational constant G = 1.) We 
use the parallel N-body tree-code GADGET (Springel, Yoshida & White 2001) to run the 
simulations. The stars in the cluster arc represented by = 10^ — 10^ equal mass particles. 
The velocities of the particles have a Maxwellian distribution. The gravitational potential 
is softened with a softening length of e = 0.77i?hA^~^^^, where i?h is the initial half-mass 
radius of the cluster. (Gravity is softened to minimize numerical noise due to two-body 
interactions between the particles.) The individual timesteps are equal to ■\/2r]e/a, where 
a is the acceleration of a particle, and the parameter r] is made small enough to ensure the 
total energy conservation to better than 1%. This configuration is the simplest realization 
of possible initial conditions for GCs, and yet results in remarkably good agreement with 
current cluster parameters. 

It is convenient to express the mass of a cluster M in units of the virial mass Myij. : M = 
M^irlO^, where (3 is a mass parameter. For a homogeneous sphere Myir = a"i^j(47rpini/375) 
The total energy of the system becomes positive (and the system becomes formally unbound) 
for p < -0.45. 
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We ran a set of 17 different models (see Table 1) with the same values of pini and 
(7ini, and parameter f3 ranging from -0.8 to 2. (The corresponding initial virial parameter 
u = 2K/W = 10~^^/^ is u = 3.4 to 0.046, where K and W are initial kinetic and potential 
energy.) Models with (3 = -0.8 are completely unbound throughout the simulation, whereas 
models with f3 = (-0.7, -0.6, -0.5) form a bound cluster, containing less than 100% of the total 
mass, after the initial expansion phase. Models with (3 > initially collapse to a smaller half- 
mass radius rmin (see Table 1), and then bounce. All models with /3 > —0.7 experience phase 
mixing and/or violent relaxation, and within 10 - 100 crossing times reach an equilibrium 
state with a flat core and sharply declining density at large radii. The projected surface 
density profiles of equilibrium clusters closely resemble the profiles of Galactic GCs. The 
match is the best for the GCs which have experienced little dynamical evolution (see Figure 

1)- 

Our models fit the observed surface density profiles of GCs very well, despite the lack 
of a tidal radius as present in King models. We do not include an external tidal field in 
our simulations, so this lack is expected. The tidal cutoff in surface density is observed in 
very few GCs (Trager et al. 1995), due to the contamination with background stars at the 
outskirts of the clusters. 

As we will show in Section 3, all real GCs should have experienced an adiabatic collapse 
(when the orbital angular momentum is conserved for individual stars). To make sure that 
our models are in the same collapse regime as the real clusters, all our models should satisfy 
the following adiabaticity criterion (Aarseth, Lin, & Papaloizou 1988): 

In our case, this condition can be expressed as /? < (1/2) log iV — (3/2)log5. The adiabatic 
collapse criterion is then /3 < 1.5 for N — 10^ and /5 < 2.0 for N — 10^. According to this 
criterion, all our models undergo an adiabatic collapse (see Table 1). In our models, the 
collapse factor C = -Rh/'^min correlates well with the virial parameter u: C oc u"^. The value 
of the exponent 7 = 0.95 ± 0.02 is very close to the adiabatic value of 7 = 1 (Aarseth et al. 
1988). 

Prom our set of models, we obtained the following correlations between the projected 
half-mass radius rh, projected central velocity dispersion ao, central surface density Eq, 
central density Pc, King core radius Tq, binding energy Ey^, and mass M (the uncertainties 
in the exponents are < 0.01): 
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Th ~ constant, ctq oc M^-^^ Eq oc M^-'^^ 

PcOcM^•6^ Eb(xM^■'^^ (2) 



These correlations are valid for (3 > 0, with the last relation being valid for (3 > —0.6. 
Most of correlations become non-linear in log-log space for /5 < (corresponding to u > 1). 
Although for our coldest models, rh scales (as expected) as the initial radius (and hence as 
M^/^), the total variation in rh over two orders of magnitude in system mass (/3 = — 2) is 
only ~ 0.12 dex (see Table 1; please note that i?ini oc M^/^). 

Another very interesting correlation is ro — rmin within the measurement errors for 
(3 > 0.2 (Table 1). The above correlation can be understood if we rewrite the theoretical 
result C = Rh/rmin oc (Aarseth et al. 1988) for the case of constant pini and fTini: 
rmin oc M~^/^. The theoretical relation is very close to the model relation ro oc M"^-'^'^, 
resulting in ro oc rmin as observed. 

We tested the numerical convergence of our results by making two additional runs for 
P — 1.6. In the first one, the number of particles was reduced to A^" = 10^, with the optimal 
value for the softening length of log(e/i?ini) = —1.88. Despite the fact that this run had the 
most poorly resolved "crunch" among all our models (rmin/e — 3.6), all the derived model 
parameters were identical to the original run parameters within measurement errors. In the 
second run, we again used a reduced number of particles A^ = 10^, but the parameter e was 
made much smaller than the optimal value: log(e/i?ini) = —2.67. In this run, the crunch 
is resolved very well — Tmin/e — 23, but again we did not see significant deviations from 
the original run parameters. The results of the additional runs suggest that the number 
of particles N and the value of the softening length e we use in our runs are adequate for 
resolving the violent relaxation process. 



3. COMPARISON WITH OBSERVATIONS 

We assume that the mass-to-light ratio M/L is a constant, which is well-established 
for observed GCs (log[M/L] = 0.16 ± 0.03, McLaughlin 2000). Then we derive the same 
correlations for observed clusters as in equations (2) using the onhne version of the Milky 
Way GC catalogue of W. E. Harris (Harris 1996). We used data for 109 non-core-collapsed 
GCs. Forty- five of these clusters have a measured ao. Binding energies were calculated using 
equations (5a-c) of McLaughlin (2000). The observed correlations are: 
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rt ~ constant, ao oc MO-^^io.os^ ^ ^i.3i±o.io^ 

Pe OC M^-53±0.17^ ^ ^-0.28±0.07^ Efe OC M^-O^iO.OS^ (3) 

The theoretical exponents differ from the observational values by one sigma or less (with the 
exception of the ao[M] correlation, where the difference is 3 sigma). 

The closeness of the model correlations (equations [2]) to the observed correlations 
(equations [3]) can be understood if we assume that in their initial non-equilibrium config- 
uration, all Galactic GCs had the same values of stellar density pini and velocity dispersion 



To estimate the values of the universal parameters pini and cxini, one can use in princi- 
ple any two or more pairs of model/observational correlations from equations (2) and (3). 
However, one has to keep in mind that few GC parameters are well suited for comparing 
the model and observed correlations. Most Galactic GCs are in advanced stages of their dy- 
namical evolution. As GCs evolve, some of their parameters deviate from initial equilibrium 
values. This process affects mainly the central core region of a cluster where encounters 
between stars are the most frequent. Analytical and numerical calculations show that at 
some point a GC should experience a runaway core collapse due to gravothermal instability 
(Spitzer 1987). Around 20% of Galactic GCs are believed to be in a post-core-collapse state 
(Harris 1996). The analytical theory of core collapse (Spitzer 1987) can be used to find 
GC parameters which are least affected by dynamical evolution. We obtained the following 
relations: 



Pe OC r^'-'\ Eo oc r^'-'\ oc r^^-'^ . (4) 

To obtain an analogous relation for binding energy E\^, we used equations (5a-c) of McLaugh- 
lin (2000). Assuming that during core collapse the tidal radius rt does not change, we derived 



^^b oc ro-°•°^ (5) 

which is accurate for the range of rt/ro = 10 — 100 corresponding to 75% of Galactic GCs. 
Equations (4) and (5) show that the parameters £'b and gq are least affected by dynamical 
evolution. Numerical simulations of core collapse show that the half-mass radius rh is also a 
very slowly evolving parameter (Spitzer 1987). 
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We used two following correlations to estimate the values of pini and ami'. £'b(co) and 
Th = constant. The fitting gave the following estimates: logpini = 1-14 ± 0.26, and 
logcTini = 0.28 ± 0.11. (The units for pini and cTini are Mq pc~^ and km s~^) The relation 
between the masses of real GCs and the model mass parameter j3 is then M ~ 3.5 x 10'^+^ Mq. 
The corresponding minimum initial cluster mass resulting in a bound cluster (the model 
with (3 = -0.7) is M^m ~ 6, 900 Mq , with a one-sigma interval of 3,000 - 16,000 Mq. The 
criterion of an adiabatic collapse of Aarseth et al. (1988) (equation [1]) can be reexpressed 
as M < 3(Tfjji/(47rpinimi), where mi is a typical stellar mass in the cluster. For our values of 
Pini and (Tini, a collapse is adiabatic if M < 1.6 x 10^ M© (we assumed mi = 0.6 Mq). As 
one can see, our adiabatic collapse models are adequate for the whole range of GC masses. 

The simplest physical interpretation of the universal initial density value is to assume 
that the major burst of star formation in a contracting proto-GC molecular cloud occurs 
when the gas density reaches a certain critical value pcr > Pm (the non-equality sign is to 
account for less than a 100% efficient star formation and mass loss due to stellar evolution). 
The initial velocity dispersion a-mi can be assumed to be commensurable with the sound 
speed in the turbulent proto-GC cloud, and hence with its temperature T, at the moment 
when the star burst occurs. For a purely molecular gas of primordial composition (Yp = 
0.291), we obtained the following estimates of the universal gas density Pcr and temperature 
T for star-forming gas in proto-GCs: pcr > 230 cm~^, and T ~ 1,000 K. The corresponding 
gas pressure isP>2.3xlO^K cm~^. 

In Figures 2-4 wc compare three different observed correlations for Galactic GCs with 
the model predictions rescaled to logpini = 1.14 and logcxini = 0.28: ii^b(c"o), -£'b(^) and 
Eo(pc). In Figures 2 and 3 (showing the correlations between the least evolved cluster 
parameters) the agreement between the model and data is excellent, with a low statistical 
significance suggestion that the most evolved clusters, shown as open circles, tend to deviate 
from the model E-^{M) relation (Figure 3). Figure 4 shows a correlation between the most 
evolved cluster parameters — central surface density Eq and central density pc (see equations 
[4]). The situation with Figure 4 is in accord with our model predictions: all clusters deviate 
from the model, "zero- age", relation in a systematic fashion, with more evolved clusters 
located farther from the model line. (Please note that the correlations from Figures 3 and 4 
were not used to fit the model to the data.) 

4. SUMMARY 

We have shown that a simple model of violent relaxation of a constant density, constant 
velocity dispersion stellar sphere produces clusters which are remarkably similar to the ob- 
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served Galactic GCs. They have the correct surface density profiles, and a large number of 
correlations between cluster parameters are accurately reproduced by our models as long as 
we assume a constant initial density and initial velocity dispersion. This result implies that 
all GCs were formed from gas with a universal value of density and temperature. In reality, 
the constant density, constant velocity dispersion initial setup of our models may correspond 
to isothermal turbulent cores of giant molecular clouds with a Gaussian density profile. Our 
simple picture of GC formation suggests that the conditions in the cluster-forming gas must 
be quite uniform across a large portion of the early universe. 

We thank W. E. Harris, R. E. Pudritz, and D. E. McLaughlin for useful discussions, 
and the anonymous referee for improvements to the paper. S. M. is partially supported by 
SHARCNet. The simulations reported in this paper were performed at CITA. 
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Table 1: Model parameters 



Input parameters 



Derived parameters 



p 


AT 
iV 


v 


Across 




1 /-V rv ^min 

R 

-'I'lni 


1„„ ro 

-'*'ini 




^og — 


^og — 


0.8 


10^ 


3.4 


5000 


-1.88 












0.7 


10^ 


2.9 


1000 


-1.88 




0.50 




-1.10 


-2.9 


0.6 


10^ 


2.5 


1000 


-1.88 




0.21 


1.09 


-0.58 


-1.3 


0.5 


10^ 


2.2 


1000 


-1.88 




0.05 


0.34 


-0.37 


-0.7 


0.4 


10^ 


1.8 


1000 


-1.88 




-0.06 


0.08 


-0.25 


-0.3 


0.3 


10^ 


1.6 


1000 


-1.88 




-0.12 


-0.06 


-0.15 


0.0 


0.2 


10^ 


1.4 


1000 


-1.88 




-0.20 


-0.15 


-0.07 


0.2 


0.0 


10^ 


1.0 


200 


-1.88 


-0.25 


-0.31 


-0.29 


0.07 


0.6 


0.2 


10^ 


0.74 


200 


-1.88 


-0.42 


-0.42 


-0.40 


0.19 


0.9 


0.4 


10^ 


0.54 


200 


-1.88 


-0.58 


-0.57 


-0.49 


0.30 


1.3 


0.6 


10'^ 


0.40 


200 


-1.88 


-0.71 


-0.70 


-0.57 


0.41 


1.6 


0.8 


10'^' 


0.29 


50 


-1.88 


-0.84 


-0.79 


-0.64 


0.52 


1.9 


1.0 


10^ 


0.22 


50 


-1.88 


-0.96 


-0.92 


-0.72 


0.63 


2.3 


1.2 


10^ 


0.16 


50 


-1.88 


-1.08 


-1.03 


-0.79 


0.74 


2.6 


1.4 


10^ 


0.12 


50 


-1.88 


-1.19 


-1.16 


-0.84 


0.85 


2.9 


1.6 


10^ 


0.086 


15 


-2.21 


-1.32 


-1.28 


-0.88 


0.98 


3.3 


2.0 


10^ 


0.046 


10 


-2.21 


-1.54 


-1.53 


-0.94 


1.20 


4.0 



Note. — Here ro, (To, and are the core parameters for relaxed models: King core radius, projected 
central velocity dispersion, and central density, respectively; rmin is the minimum half-mass radius during 
the collapse, and rh is the projected half-mass radius for relaxed models. 
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Fig. 1. — Surface density profile (solid line) for a model with P — 0.4. This model matches 
well the surface density profiles of five GCs (points, Tragcr, King & Djorgovski 1995) which 
are among the least dynamically evolved ones. Surface density is normalized to 1 at the 
center of the cluster. The vertical short-dashed and long-dashed lines show the values of the 
softening length e and the model core radius ro, respectively. 




Fig. 2. — Cluster binding energy as a function of velocity dispersion for 45 Galactic GCs with 
a measured o"o (circles) and our model (solid line). The open circles are the most dynamically 
evolved clusters with logtc < 8.3, where tc is the core relaxation time in years. 
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Fig. 3. — Cluster binding energy as a function of total cluster mass for 109 non-core-collapsed 
GCs (circles) and our model (solid line). The open circles are those clusters with logic < 8.3 
(as in Figure 2). 



- 12 - 




Fig. 4. — Central surface density as a function of central density for 109 non-core-collapsed 
GCs (circles) and our model (solid line). The open circles are more dynamically evolved 
clusters with logic < 9-2. 



